05. stat & ML basic

Author

GC

Published

May 12, 2024

모집단과 표본

1 모집단(Population) : 연구 또는 관측 대상이 되는 전체

  • ex : 우리나라 30세이상 월 평균 소득 조사, 사실 이를 위해 해당 집단을 전수 조사하는 것은 거의 불가능에 가까움

2 표본 : 모집단의 일부, 연구 또는 관측을 위해 일정 시점 또는 일정 집단에서 수집한 대상(데이터)

  • ex : 2020년 1월 1일시점에 조사한 30세 이상 국내 성인 남자 10명의 월평균 소득

3 표본평균과 중앙값 구하기

  • 평균
income <- c(100, 200, 300, 400, 500, 600, 700, 800, 1000 ,1000)

sum(income)/length(income)
560
  • 중앙값
#sort(income);length(income)

print(paste("손계산 :",(sort(income)[5]+sort(income)[6])/2))

print(paste("함수 사용 :",median(income)))
[1] "손계산 : 550"
[1] "함수 사용 : 550"

평균과 분산

- 모집단의 분포(형태는) 일반적으로 수집한 표본의 모수인 평균과 분산에 의해 결정된다.

평균(기댓값)

1 이산형일 경우 : 베르누이, 이항, 포아송 분포 등등…

ex1 : 동전 던지기를 해서 앞면이 나올 평균확률은?

\(x\) 0 1
\(f(x)\) 1/2 1/2

\[E(X) = \sum_{i=1}^{n} x_i f(x_i)\]

2 연속형일 경우 : 베타, 지수, 감마, 정규 등등…

ex2, 표준정규분포의 평균은?

\[X \sim N(0, 1)\]

\[f(x) = \frac{1}{\sqrt {2\pi}} e^{- \frac {x^2}{2}}\]

숙제 1. 아래의 공식을 이용하여 표준정규분포의 평균이 실제로 1이 나오는지 확인!

\[E(X) = \int xf(x) \,\,dx\]

대수의 법칙

- Law of large numbers

- 샘플사이즈가 클수록 표본평균은 모집단의 평균으로 수렴한다는 법칙! (체비셰프 부등식을 이용하여 증명 가능!)

\[\begin {align} P(|\overline X - \mu| < \varepsilon) &= P(|\overline X - \mu| < \varepsilon^2) \\ \\ &= 1 - \frac{E(\overline X_{n} - \mu)^2}{\varepsilon^2} \\ \\ &= 1 - \frac{\sigma^2/n}{\varepsilon^2} \approx 1 \end {align} \]

- R을 이용한 증명

\(x\) 100 200
\(f(x)\) 1/2 1/2
  • 이므로 모집단의 기대값은 $E(X) = xf(x) = 100 /2 + 200 /2 = 150 $
### code_fold : true
#install.packages("plotly")
library(tidyverse)
library(plotly)


x <- tibble(income = c(100, 200))

data <- tibble()
for (i in 1:500) 
    {
        sample <- x %>% sample_n(i, replace = T, weight = c(1/2,1/2))
        sample <- sample %>% mutate(sample_size = i)
        data <- bind_rows(data, sample)
    }
### code_fold : true
sample_mean <- data %>% group_by(sample_size) %>% 
                        summarize(sample_mean = mean(income))

fig <- plot_ly(data = sample_mean, 
               x = ~sample_size, y = ~sample_mean, mode = "lines", 
                type = "scatter", name = "표본평균", opacity = 0.5) %>% 
                        add_trace(y = 150, name = "모집단 평균")

fig %>% 
    layout(title = "Law of large numbers") 

분산

- 수집한 표본들이 평균으로부터 떨어져 있는 정도를 측정

1 이산형일 경우

ex1 : 동전 던지기를 해서 앞면이 나오는 실험의 분산은?

\(x\) 0 1
\(f(x)\) 1/2 1/2

\[\sum_{i=0}^{n} = [x_i-E(X)]^2f(x_i)\]

2 연속형일 경우

\[Var(X) = \int [x-E(x)]^2f(x) \,\,dx\]


정규분포

1 분포 정의 \(\to X \sim N(\mu, \sigma^2)\)

\[ f(x) = \frac{1}{\sqrt(2\pi\sigma^2)}\exp\left( -\frac12 \left(\frac{x-\mu}{\sigma}\right)^2\right)\]

2 대게 종모양의 분포를 띔

3 표준편차가 작을수록 곡선이 좀더 뾰족해지고 폭이 좁아진다? \(\to\) 분산에 정의를 다시 한번 생각해보시면됩니다.

options(repr.plot.res=200, repr.plot.width = 8, repr.plot.height = 6)
x1 <- rnorm(10000, 0, 1)
x2 <- rnorm(10000, 0, 10)
x3 <- rnorm(10000, 0, 20)
x4 <- rnorm(10000, 0, 30)

par(mfrow = c(2,2))
hist(x1, main = "평균 : 0, 표준편차 : 1", breaks = 50, xlim = c(-100,100))
hist(x2, main = "평균 : 0, 표준편차 : 10", breaks = 50, xlim = c(-100,100))
hist(x3, main = "평균 : 0, 표준편차 : 20", breaks = 50, xlim = c(-100,100))
hist(x4, main = "평균 : 0, 표준편차 : 30", breaks = 50, xlim = c(-100,100))

4 정규분포의 기각영역

  • 우리가 수집한 어떤 집단의 평균이 0일 때의 가설검정

\[H_0 : \mu = 0, \quad H_1 : \mu \neq 0\]

x = rnorm(1000)
normal <- tibble(x = x, fx = fx)

normal <- normal %>% mutate(p = ifelse(abs(x) >= 2.58, "유의수준 1%", 
                             ifelse(abs(x) >= 1.96, "유의수준 5%", 
                                    ifelse(abs(x) >= 1.65, "유의수준 10%", "기각 x"))))
options(repr.plot.res = 200, repr.plot.width = 8, repr.plot.height = 4)
normal %>% ggplot(`aes(x = x, fill = p)) +
            geom_histogram(bins = 100)

5 중심극한정리

  • 표본의 크기가 커지면, 표본평균의 분포는 정규분포로 수렴 \(n \to \infty,\, \overline X_n \sim N(\mu, \sigma^2/n)\)

  • 이것을 가정하고 있기 때문에, 우리는 어떠한 가설설정 시 정규분포를 이용해서 통계적 유의성을 말할 수 있는거임

example1, 동전던지기를 통한 중심극한정리 실험

  • step1. 동전던지기 실험, 동전던지기 100회 \(\text{sample size} = 100\)
rbinom(100, 1, 0.5)
  1. 0
  2. 1
  3. 1
  4. 0
  5. 1
  6. 1
  7. 0
  8. 1
  9. 1
  10. 1
  11. 1
  12. 1
  13. 1
  14. 1
  15. 1
  16. 0
  17. 0
  18. 1
  19. 1
  20. 1
  21. 0
  22. 1
  23. 1
  24. 0
  25. 0
  26. 0
  27. 0
  28. 0
  29. 1
  30. 0
  31. 0
  32. 1
  33. 0
  34. 0
  35. 0
  36. 1
  37. 1
  38. 1
  39. 0
  40. 0
  41. 1
  42. 0
  43. 1
  44. 1
  45. 1
  46. 1
  47. 1
  48. 0
  49. 1
  50. 0
  51. 0
  52. 0
  53. 0
  54. 0
  55. 0
  56. 1
  57. 1
  58. 0
  59. 1
  60. 1
  61. 1
  62. 1
  63. 0
  64. 0
  65. 0
  66. 0
  67. 0
  68. 1
  69. 0
  70. 1
  71. 0
  72. 0
  73. 0
  74. 0
  75. 1
  76. 1
  77. 1
  78. 0
  79. 0
  80. 0
  81. 0
  82. 1
  83. 0
  84. 1
  85. 1
  86. 1
  87. 1
  88. 1
  89. 0
  90. 1
  91. 0
  92. 1
  93. 1
  94. 1
  95. 0
  96. 0
  97. 0
  98. 0
  99. 0
  100. 0
  • step2. 표본평균 구하기
x_mean = mean(rbinom(100, 1, 0.5))
x_mean
0.6
  • step3. 표본평균의 분포 구하기
sample_mean = c()

for (i in 1:1000) {
    x_mean = mean(rbinom(1000, 1, 0.5))
    sample_mean = append(sample_mean,x_mean)
    }
hist(sample_mean,main = "n=1000", breaks = 50)


ML basic

1. 선형회귀 평가지표

- 결정계수 \(R^2\) : 우리가 적합시킨 모델이 얼마만큼의 설명력을 가지고 있는가

  • 지난시간 학습한 자료 꼭 한번 확인!

\[R^2 = \frac{SSR}{SST} = \frac{1-SSE}{SST}\]

2. Overfitting

  • 우리가 학습시킨 모델이, 훈련데이터에만 집중적으로 학습이되어서, 실제 검증용 데이터로 평가를 했을 때 모델 성능이 낮아지는 경우

  • 이거를 표현해보면…

\[\begin {align} \text {MSE} = E[(\hat y - y)^2] &= Var(\hat y) + \text {bias}^2 \\ &= E[(\hat y - E(\hat y ))^2] + E[(E(\hat y) - y)^2] + Var(\varepsilon) \end {align}\]

- \(Var(\hat y) \, \to \,\) 훈련자료의 변화에 \(\hat y\)가 얼마나 민감하게 반응하는가?(즉, 설명력을 의미한다.)

- \(\text {bias}^2\) : 실제자료를 모형으로 얼마나 가깝게 근사할 수 있는가? (만약, \(E(\hat y) = y\) 일 경우 해당 추정량은 불편추정량이다.)

  • \(Var(\hat y),\,\text {bias}^2\) 는 reducible error로 어떤 모델을 선택하느냐에 따라 줄일 수 있는 오차이다.

- \(Var(\varepsilon)\) : 어떤 모델을 사용하건 줄일 수 없는 오차

- Flexibility : 모델 복잡도 측도

  • Flexibility가 높아지면 \(Var(\hat y)\) 는 높아지고, \(\text {bias}\)는 낮아진다. (trade-off)

  • 이를 바꿔말하면 더 유연한(복잡한) 모형은 더 큰 분산을 가지고 더 단순한 모형일수록 큰 편의를 가진다.

  • 훈련 데이터로 모델 적합시, bias를 낮추기 위해 Var를 높여 모델 복잡도를 지나치게 높이는 것은 좋은 선택이 아니다.

  • 아래의 예시는 모형의 복잡도가 너무 높으면 평가 데이터에 대한 MSE가 오히려 과대추정되는 “과적합 문제” 를 보여주고 있다.

ex1

- 실제로 위의 원데이터는 일반적인 선형모형으로 적합해도 괜찮을 것 같다.

- 왼쪽그림에서 지나치게 복잡하게 적합된 초록색 선을 살펴보자.

  • 훈련 MSE(회색선) 는 3개의 모델 중 가장 낮으나, 과적합 이슈로 평가 MSE(빨간색선) 가 과대추정된 것을 볼 수 있다.

ex2

- 위 데이터는 ex1과 달리 일반적인 선형회귀모형으로 적합하면 안될것 같은 느낌이 든다.

- 파란색선은 중간 정도로 적합된 복잡한 모형이고, 초록색선은 지나치게 적합된 모형이다.

  • 오른쪽 그림을 살펴보면 파란색선으로 적합된 모형의 Test MSE가 가장 낮음을 볼 수 있다.

과적합 방지 1. 변수선택법

- 변수선택법 : AIC, BIC, Mallows’ C_p 등등

과적합 방지 2. Regularized Linear Regression

- 정규화 선형회귀는 선형회귀 계수에 대한 제약 조건을 추가하여 모델이 과도하게 최적화되는 현상을 막는 방법!

1 Ridge

  • 기존 OLS를 이용한 \(\beta_i\) 추정할 떄 사용하는 정규방정식

\[RSS = \sum_{i=1}^n(y_i - \beta_0 -\sum_{j=1}^p\beta_jx_{ij})^2\]

  • Ridge 추정량은 위와 매우 유사하나 shrinkage penalty를 추가로 고려하여 최소화한다.

\[\sum_{i=1}^n(y_i - \beta_0 -\sum_{j=1}^p\beta_jx_{ij})^2 + \lambda \sum_{i=1}^p \beta_j^2\]

\[\lambda \,:\,tuning \,\,parameter\]

  • \(\lambda\)가 0이되면 일반적인 선형회귀모형이되고, \(\lambda\)가 커지면 정규화 정도가 커진다.(회귀계수들이 작아진다.) \(\to\) bias는 증가 variance는 감소

  • 릿지회귀모형에서는 L2 penalty 를 사용하여 베타계수들이 0에 근사하도록 한다.

2 Lasso

  • 라쏘는 릿지와 다르게 베타계수를 “0값으로 보낸다.”

  • 이것을 “L1 penalty” 라고 부른다.m

\[\sum_{i=1}^n(y_i - \beta_0 -\sum_{j=1}^p\beta_jx_{ij})^2 + \lambda \sum_{i=1}^p |\beta_j|\]

실습

#install.packages("mosaicData")
library(tidyverse)
library(mosaicData)
library(glmnet) ## 리지 및 라쏘 회귀를 적합하기 위한 패키지

1 데이터확인

glimpse(RailTrail)
Rows: 90
Columns: 11
$ hightemp   <int> 83, 73, 74, 95, 44, 69, 66, 66, 80, 79, 78, 65, 41, 59, 50,~
$ lowtemp    <int> 50, 49, 52, 61, 52, 54, 39, 38, 55, 45, 55, 48, 49, 35, 35,~
$ avgtemp    <dbl> 66.5, 61.0, 63.0, 78.0, 48.0, 61.5, 52.5, 52.0, 67.5, 62.0,~
$ spring     <int> 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0,~
$ summer     <int> 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,~
$ fall       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1,~
$ cloudcover <dbl> 7.6, 6.3, 7.5, 2.6, 10.0, 6.6, 2.4, 0.0, 3.8, 4.1, 8.5, 7.2~
$ precip     <dbl> 0.00, 0.29, 0.32, 0.00, 0.14, 0.02, 0.00, 0.00, 0.00, 0.00,~
$ volume     <int> 501, 419, 397, 385, 200, 375, 417, 629, 533, 547, 432, 418,~
$ weekday    <lgl> TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TR~
$ dayType    <chr> "weekday", "weekday", "weekday", "weekend", "weekday", "wee~

- model.matrix는 기본적으로 절편을 포함한 모형을 산정한다. 절편을 제외시킨 모형을 고려해보자

x <- model.matrix(volume~.-1 ,RailTrail) 
y <- RailTrail$volume

- 모형적합

ridge.fit <- cv.glmnet(x,y,alpha=0) ## 모형적합 alpha=0 은 릿지를 말함
bestlam <- ridge.fit$lambda.min
log(bestlam)
3.02271971457234

- 각 변수들의 beta 계수구하기

coef(ridge.fit, s = "lambda.min")
12 x 1 sparse Matrix of class "dgCMatrix"
                        s1
(Intercept)    111.4030636
hightemp         3.6351375
lowtemp         -0.9192012
avgtemp          1.9016337
spring           9.6677301
summer           4.4014186
fall           -27.7525473
cloudcover      -8.1350693
precip         -88.4982099
weekdayFALSE    13.8267092
weekdayTRUE    -13.8900786
dayTypeweekend  13.9990402

- MSE가 최소가 되는 지점을 확인

options(repr.plot.res=200,repr.plot.height=4,repr.plot.width=10)

plot(ridge.fit)
abline(v=log(bestlam),lty="dashed",lwd=2,col="blue")

- 산출된 best.lam값으로 모형을 다시 적합후 test_MSE를 구해보자

train  <- sample(1:nrow(x),nrow(x)*0.7)
test  <- -train
y.test <- y[test]
y.train <- y[train]
ridge.fit <- glmnet(x[train,], y.train, alpha=0, lambda=bestlam, family="gaussian")

ridge.pred <- predict(ridge.fit,s=bestlam,newx=x[test,]) 

ridge.coef <- predict(ridge.fit,s=bestlam,newx=x[test,],type="coefficients") ## type="coefficients"로 하면 예측된 베타계수를 보여줌

- 여기선 MSE가 거의 1만에 가깝게 나왔는데 이는 지나치게 모형을 단순화한 것임 즉, best_lam값을 다시 찾는 과정을 반복해야한다….

mean((ridge.pred - y.test)^2)
8763.80834274373